Predicting the current fishable habitat distribution of Antarctic toothfish (Dissostichus mawsoni) and its shift in the future under climate change in the Southern Ocean

Global warming continues to exert unprecedented impacts on marine habitats. Species distribution models (SDMs) are proven powerful in predicting habitat distribution for marine demersal species under climate change impacts. The Antarctic toothfish, Dissostichus mawsoni (Norman 1937), an ecologically and commercially significant species, is endemic to the Southern Ocean. Utilizing occurrence records and environmental data, we developed an ensemble model that integrates various modelling techniques. This model characterizes species-environment relationships and predicts current and future fishable habitats of D. mawsoni under four climate change scenarios. Ice thickness, depth and mean water temperature were the top three important factors in affecting the distribution of D. mawsoni. The ensemble prediction suggests an overall expansion of fishable habitats, potentially due to the limited occurrence records from fishery-dependent surveys. Future projections indicate varying degrees of fishable habitat loss in large areas of the Amery Ice Shelf’s eastern and western portions. Suitable fishable habitats, including the spawning grounds in the seamounts around the northern Ross Sea and the coastal waters of the Bellingshausen Sea and Amundsen Sea, were persistent under present and future environmental conditions, highlighting the importance to protect these climate refugia from anthropogenic disturbance. Though data deficiency existed in this study, our predictions can provide valuable information for designing climate-adaptive development and conservation strategies in maintaining the sustainability of this species.


INTRODUCTION
Anthropogenic-induced climate changes are continuing to alter the marine ecological environment, such as an increase in water temperature, ocean acidification and a decrease of ice cover in the polar ocean (Perovich & Menge, 2009;Duncan et al., 2022;Falcón & niche of species (Citores et al., 2020), estimating the risk of invasive species (Carman et al., 2016;Poursanidis et al., 2020), and informing conservation strategies (Yoon et al., 2017;Cázares et al., 2021) over the past two decades.The established SDMs based on the current data can forecast the future distribution of species on condition that future environment conditions can be projected based on the global climate models under different greenhouse gas emission scenarios.SDMs are constructed using various statistical algorithms.Integrating these into a multi-algorithm ensemble model overcomes the limitations of individual models, improving predictive accuracy and reducing uncertainty (Robinson et al., 2017).Recently, a growing number of studies have applied the ensemble model to estimate climate impacts on the habitat suitability of demersal deep-sea species, such as Antarctic jonasfish and Pacific cod (Li et al., 2022;Ran et al., 2022).Yates et al. (2019) integrated fishery-dependent and environmental data to investigate the distribution of D. mawsoni along East Antarctica.However, how the climate change to affect the habitat suitability of D. mawsoni under the mid-and long-term effects did not attract enough attention at present.
In this study, we first obtained occurrence records and environment data from online public datasets; we then developed an ensemble model by integrating individual models to characterize species-environment relationships; we finally conducted the ensemble predictions of D. mawsoni fishable distribution under present-day and future environment conditions in the Southern Ocean.This study aims to: (1) identify key factors influencing D. mawsoni's distribution; (2) assess its fishable distribution shifts by mid and late 21st century under four climate scenarios; (3) identify climate refugia by comparing the fishable habitat distributions of D. mawsoni under current and future environmental conditions.We anticipate our study will inform climate-adaptive conservation and management strategies to sustain D. mawsoni.

Presence/pseudo-absence data
The occurrence records of D. mawsoni were downloaded from three public online datasets: the Ocean Biodiversity Information System (OBIS, http://www.iobis.org), the Global Biodiversity Information Facility (GBIF, https://doi.org/10.15468/dl.h33snm),and the Integrated Digitized Biocollection (iDigBio, http://www.idigbio.org).A total of 11,849 occurrence records located in the Southern Ocean between 2000-2014, which corresponded with the time span of environment variables, were obtained.These records were mainly from fishery-dependent sampling and therefore our predictions based on them can only represent the fishable distribution of D. mawsoni.Because sampling was not random, the spatially clustering of occurrence record may over-represent the environment conditions.To mitigate this bias, we used the R package spThin (Lammens et al., 2015) to filter records, retaining one per 5 × 5 arc min grid (approximately 9.2 × 9.2 km at the equator), aligning with the environmental data's spatial resolution.After this data thinning process, 451 occurrence records were retained for the model analysis (Fig. 1).
Following Iturbide, Bedia & Gutiérrez (2018), we used a three-step method to generate pseudo-absence records, tenfold the number of occurrence records.The initial step involved utilizing a presence-only support vector machine algorithm to define environmentally unsuitable regions and then limited pseudo-absence records to these unsuitable regions.The subsequent step consisted of constructing SDMs using pseudo-absences within buffers of varying extents encircling the occurrence sites.The third step focused on identifying the optimal buffer distance informed by the performance of the SDMs established in the preceding step.This method balances the use of spatial and environmental dimensions for the selection of pseudo-absence points and has been proven to be superior to other methods.

Environmental predictor variables
Selecting appropriate predictor variables is crucial for the effectiveness of Species Distribution Models (SDMs).Considering data availability and species-environment relationships, ten environmental predictors were initially selected.Mesa, Riginella & Jones (2019) found that juvenile D. mawsoni was not uniformly distributed along the East Antarctica and its spatial variations involved complex interactions with bathymetry, temperature, sea ice and salinity.These variables were thus chosen as predictors.As a high-trophic level predator, the spatial distribution of D. mawsoni is greatly affected by the variations in the abundance and distribution of prey.D. mawsoni is not herbivorous so chlorophyll concentration cannot affect its distribution directly, but this factor may affect the prey availability through the bottom-up effects.Therefore, this factor was involved in the model analysis.D. mawsoni, predominantly found on continental shelves and slopes, has a circumpolar distribution.Mature adults prefer deeper waters for spawning (Hanchet et al., 2015), making seabed slope, ruggedness, and Bathymetric Position Index (BPI) relevant to their distribution.These factors were calculated from the depth layer utilizing ArcGIS Pro v3.0 and then were incorporated into our model analysis.The slope grid was generated using the integrated Spatial Analyst tool.Ruggedness and BPI were calculated using the Benthic Terrain Model 3.0 tool.For ruggedness, a neighborhood size of five grid cells was applied, while BPI was calculated with an inner radius of one and an outer radius of five grid cells, in accordance with the methodology outlined by Walbridge et al. (2018).
Ocean circulation primarily affects the early life-history stages by transporting eggs, larvae, and juvenile Antarctic toothfish from the offshore spawning areas to the coastal recruitment areas (Behrens et al., 2021).It also facilitates adults to swim along the shelf slope and back into the spawning grounds (Ashford et al., 2012).Therefore, we included current velocity as a predictor.
Multicollinearity among predictors can lead to model over-parameterization, compromising predictive accuracy (Ran et al., 2022).To address this, we conducted a variance inflation factor (VIF) analysis.The predictor manifesting a value of VIF over 10 was systematically excluded (O'Brien, 2007).After this process, five terrain variables (depth, distance to shore, slope, ruggedness and BPI) and five temporally dynamic environmental variables (current velocity, salinity, mean temperature, mean chlorophyll concentration and ice thickness) were retained for the model analysis (Table 1).Additionally, these variables exhibited low paired Pearson correlation coefficients (Fig. 2).
The Bio-ORACLE v2.2 dataset (http://www.bio-oracle.org)(Assis et al., 2018) offers environmental data for current (2000-2014 averages) and future periods (2050s: 2040-2045 averages, 2100s: 2090-2100 averages) under various greenhouse gas emission scenarios, with a 5 × 5 arc min spatial resolution.The future environment data were projected based on three atmospheric-ocean general circulation models (AOGCMs: CCSM4, HadGEM2-ES and MIROC5).In order to reduce model-dependent uncertainty, we used the average values of predictions from the three AOGCMs to represent the future environment conditions.We considered four Representative Concentration Pathways (RCPs): RCP 2.6, RCP 4.5, RCP 6.0 and RCP8.5, which indicate the continually increased  concentrations of greenhouse gas emissions.Depth and distance to shore, sourced from the Global Marine Environment Datasets (Basher, Bowden & Costello, 2018), were assumed to remain constant over time.

Modelling process
Utilizing presence/pseudo-absence records and current environment data, we developed SDMs to relate the occurrence probability of species with environment conditions using the R package biomod2 (Thuiller et al., 2014).The package incorporates ten widely-utilized modeling algorithms: Surface Range Envelope (SRE), generalized linear model (GLM), generalized additive model (GAM), multiple adaptive regression splines (MARS), generalized boosting model (GBM), flexible discriminant analysis (FDA), classification tree analysis (CTA), Random Forest (RF), artificial neural network (ANN), and maximum entropy (Maxent).The predictive ability of each model was evaluated by a five-fold cross validation method in which 80% of data were selected randomly to train the model and the rest 20% of data used to validate the model (Guisan, Thuiller & Zimmermann, 2017).This procedure was iterated 10 times.Area under the curve (AUC) and true skill statistics (TSS) were used to assess model performance.To reduce model-dependent uncertainty, we used the AUC value of each individual model as weight to build an ensemble model.Models exhibiting TSS > 0.8 and AUC > 0.9 were chosen to build the ensemble model (Araújo et al., 2005).

Predictor importance and response curves
The importance of each predictor was evaluated using the random isolation method included in the biomod2 package.This method involved calculating the Pearson correlation coefficient between predictions made using all predictors and those using a randomly permuted evaluated predictor (Guisan, Thuiller & Zimmermann, 2017).High correlation signified low importance of the evaluated predictor.Response curves were plotted to visualize how occurrence probability varied with each environmental predictor's gradient.

Predictions of current and future fishable D. mawsoni distribution
The ensemble model developed based on the current environment data was used to predict the current fishable D. mawsoni distribution and future fishable distribution under four (RCP 2.6, RCP 4.5, RCP 6.0 and RCP8.5) climate change scenarios in the 2050s and 2100s.
The prediction from the model at each raster grid indicates the habitat suitability with a value between 0-1.To facilitate interpretation, continuous habitat suitability values were converted to binary (0/1) values.The grids with the habitat suitability greater than the cut-off value which maximizes the probability threshold of TSS were converted to 1 and the others were converted to 0 (Guisan, Thuiller & Zimmermann, 2017).Using the current fishable habitat distribution as a reference, the future habitat change in each raster grid has three situations: present now and in the future (stable), present now but predicted to be lost in the future (loss), and absent now but predicted to be present in the future (gain).

Model performance
TSS values ranged from 0.51 (SRE) to 0.87 (GBM), while AUC values varied from 0.75 (SRE) to 0.98 (RF) across the ten models.Both TSS and AUC revealed significant variations in predictive ability among the models (Fig. 3).Five models-GBM, GAM, CTA, RF, and MAXNET-were chosen for the ensemble model based on their TSS and AUC cut-off values.The TSS and AUC values of the ensemble model were 0.922 and 0.994 respectively, greater than any of the individual models, indicating an improvement of model performance.

Predictor importance and response curves
Ten environment predictors were included in the ensemble model to predict the potential fishable D. mawsoni distribution.Based on a randomization method, ice thickness was proved to be the most influential factor to the D. mawsoni distribution, followed by depth.Mean temperature, distance to shore, salinity and mean chlorophyll concentration showed a similar, moderate importance and the remained factors including slope, ruggedness, BPI and current velocity presented little influence on the D. mawsoni fishable distribution (Fig. 4).Response curves depicting the relationship between occurrence probability and environmental gradients are presented in Fig. 5. Occurrence probability sharply increased with ice thickness, peaking at 0.25-0.75m, before gradually declining.Occurrence probability increased as depth decreased, stabilizing at approximately 0.73 in depths of 1,000-2,000 m.Occurrence probability surged with rising water temperature, peaking at 1.5 C, and then sharply declined, bottoming out at temperatures above 5 C. Salinity, mean chlorophyll concentration, current velocity, distance to shore, ruggedness, and BPI each differently affected the fishable distribution, whereas slope changes had no discernible impact on D. mawsoni's occurrence probability.

Present-day fishable habitat distribution
The prediction from the ensemble model under current environment condition (Table 1) showed that the fishable habitats of D. mawsoni concentrated around the peripheral seas of the Antarctic continental shelf (Fig. 1).Substantial suitable fishable habitats were located in the Amundsen Sea, Ross Sea, Amery Ice Shelf, and Bellingshausen Seas.Areas predicted to have high habitat suitability encompassed most occurrence points.Relatively low habitat suitability was observed in the coastal waters of Northern Antarctica, particularly in the Weddell Sea.

Future fishable habitat distribution and its shift under future environment conditions
Environmental conditions for the 2050s and 2100s under four emission scenarios are presented in Table 2.Under all emission scenarios, current velocity and salinity exhibited minimal changes in both the 2050s and 2100s.As expected, mean water temperature showed a substantial rise and ice thickness showed a substantial decrease in the future.The extent of these changes intensified with higher emission concentrations.Unexpectedly, chlorophyll concentration increased in future projections across all emission scenarios.The fishable habitat suitability in the 2050s and 2100s under the four emission scenarios were shown in Figs. 6 and 7.The areas with different habitat suitability were distributed in the nearshore waters encircling the Antarctic continent.In the 2050s, substantial suitable habitats were located near the Antarctic Peninsula, in the Bellingshausen Sea, the Ross Sea, and around the Amery Ice Shelf.The extent of highly suitable areas expanded in 2010, with the degree of expansion varying according to the emission scenarios.Compared with the current distribution of D. mawsoni, the ensemble model predicted a medium-term (2050s) expansion of suitable habitat ranging from +9.35% (RCP 4.5) to +12.38% (RCP 8.5), and a long-term (2100s) expansion from +13.48% (RCP 2.6) to +19.81% (RCP 6.0) (Table 3).Predictions indicate that D. mawsoni's fishable habitats will extend to new areas near the Weddell Sea and the Ross Sea (Figs. 8 and 9).The extent of this expansion varies with time and greenhouse gas emission scenarios.It is noteworthy that the range of suitable fishable habitat around the Amery Ice Shelf will decrease gradually in the future (Figs. 8 and 9).Under all emission scenarios, suitable habitats in the Ross Sea are expected to persist through the 2050s and 2100s.

DISCUSSION
Ongoing global warming continues to affect marine communities and alter ecosystem functions (IPCC, 2022;Trathan & Agnew, 2010;Wassmann et al., 2011).Fish endemic to the Southern Ocean are potentially more vulnerable to global climate change due to their narrow ecological niches and limited colonization space (Freer et al., 2019).Consequently, we aimed to assess the impact of climate change on D. mawsoni's habitat suitability by integrating species occurrence records with environmental data, to guide sustainable development and management strategies.

Model performance and limitations
The ensemble model consistently outperformed individual models, as indicated by TSS and AUC values, suggesting that integrating individual models to create a weighted average enhances predictive capabilities.Previous studies corroborate the superiority of ensemble models over individual models in evaluating climate change impacts on marine demersal species (Behrens et al., 2021;Li et al., 2022;Ran et al., 2022).In this study, the SRE model exhibited the poorest performance, likely because it solely relies on presence data, unlike other models which utilize both presence and absence data.This underscores the need for simulating pseudo-absence records to enhance model performance in the absence of true absence data.
Despite the aforementioned advantages of the ensemble model, there were limitations in the data used for its development.First, the occurrence records used in this study, sourced from online public datasets, predominantly relied on fishery-dependent observations.Therefore, our study's predictions primarily represent the fishable distribution of D. mawsoni.Additionally, the fishery-dependent data held by CCAMLR, owned by  multiple countries, necessitates obtaining permission from all data owners for specific analyses.It is hard for us to gain permission for using CCAMLR data to conduct our analysis.The limited number of occurrence records may inadequately represent D. mawsoni's realized niche, thus impacting the SDM's predictive accuracy.Second, treating simulated pseudo-absence records as true absences poses certain risks.Pseudo-absences might inaccurately represent true absences, potentially due to insufficient sampling at certain sites.When true absence data are available, the SDM predictions based on presence/absence data are more accurate than those based on presence only or presence/pseudo-absence data (Xu et al., 2022).Third, given that adult D. mawsoni predominantly inhabit the sea bottom, sediment characteristics like median grain size and type may influence their distribution.Meanwhile, as a commercial-target species, human fishing activities could also affect the distribution of D. mawsoni.Failing to include substrate and anthropogenic variables due to data unavailability in the SDM will minimize the model predictive ability.Fourth, inter-annual environmental variations, such as the timing of ice emergence and retreat, influence D. mawsoni recruitment (Behrens et al., 2021).However, due to data constraints, we used multi-year (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) average environmental data to represent current conditions, thereby obscuring inter-annual variations and potentially reducing SDM predictive accuracy.Fifth, it is essential to recognize that dispersal ability significantly affects D. mawsoni's distribution.Lenoir et al. (2020) discovered that marine species adapt to suitable niches in response to climate warming six times faster than terrestrial species, attributed to the less constrained physical environment for dispersal in marine habitats.D. mawsoni is not a powerful swimmer and actually, most individuals move few kilometers throughout their entire lives (CCAMLR Secretariat, 2023).Variable ocean currents and complex bottom topography may impede this species from reaching all suitable habitats.Consequently, our study might overestimate D. mawsoni's potential future fishable distribution under the assumption of unlimited dispersal.

Response of D. mawsoni to environment change
The species-environment relationships characterized in this study can provide key information to understand the response of D. mawsoni to environment change.
Ice thickness is the most important factor in shaping the distribution of D. mawsoni.Parker et al. (2021) indicated the significant role of sea ice in D. mawsoni's early life history.Consequently, global warming-induced changes in ice thickness will substantially impact this species' population size and distribution.Our findings reveal that D. mawsoni's occurrence probability sharply increases with ice thickness values below 0.25 m and peaks at 0.25-0.75m, suggesting that specific sea ice thickness levels are vital for the species.
While the upper limit of optimal ice thickness for D. mawsoni remains unconfirmed, the reduced occurrence probability above 0.75 m could be attributed to sampling bias, as thicker ice impedes fishing boats' access, resulting in fewer captures.Yates et al. (2019) found a similar shape of response curve with ours in studying the relationship between catch rate of D. mawsoni and mean ice cover based on fishery-dependent and environment data from years 2003-2017 along East Antarctica.Ice thickness is most likely to positively correlate with mean ice cover and the similar findings confirmed the importance of sea ice for D. mawsoni.This study proved that water depth was the second important factor in shaping the distribution patterns of D. mawsoni.The previous studies in the Ross Sea and the eastern Antarctica demonstrated that the mean weight and proportion of adults increased with the increasing depth, implying a gradual migration from shallow to deep waters as fish grows (Mesa, Riginella & Jones, 2019).The Ross Sea contained important spawning grounds and the older and larger fishes spawn on the ridges, banks and seamounts of the Pacific-Antarctic Ridge to the north and east of the Ross Sea with deeper waters (Parker et al., 2019).The southern part of the Amundsen Sea was considered to be an important source of immature fishes based on the analysis of fishery data in this area (Abecasis, Afonso & Erzini, 2014;Mekonnen & Hoekstra, 2014).The acquired occurrence records were mainly located in the west of Antarctic Peninsula, the Amundsen Sea and the Ross Sea, resulting in a larger proportion of juvenile records.The preference of younger cohorts for shallower waters likely explains the observed increase in occurrence probability with decreasing depth in this study.Future work should focus on characterizing species-environment relationships, taking into account sex and ontogenetic stages.
Habitat change in the future and implication for management Beyond our expectation, a fishable habitat expansion was predicted under all emission scenarios in the future.We speculate this result may stem from the limited scope of available occurrence records, primarily derived from fishery-dependent surveys, covering only a fraction of D. mawsoni's distribution.The occurrence records were primarily located in the coast of West Antarctica, while fishing activities were also conducted in the East Antarctica in addition to in the Ross and Amundsen Sea of the West Antarctica.Due to the harsh weather conditions in the Southern Ocean, the extensive fishery-independent surveys were restricted by financial and logistical terms.Additionally, while regions with thicker ice may be suitable for D. mawsoni, such ice impedes fishing boats' access.The lack of sampling in the habitats with thicker ice is probably the reason for the predictive results of future habitat expansion.Although the data deficiency, our predictions can provide valuable knowledge for informing conservation measures and management strategies in maintaining the sustainability of D. mawsoni.The climate of the Antarctic Peninsula is the most rapidly changing in the Southern Hemisphere and marine species in this region are extremely sensitive to their environment, with predictive species and population removal in the face of a little rise in the water temperature (Mao et al., 2021).Predictions for West Antarctica were more robust due to the greater data availability in this region.Large amounts of suitable habitats in the West Antarctic Peninsula were predicted to be lost in the face of future climate change, which was consistent with a similar study on the Antarctic jonasfish (Ran et al., 2022).Parker et al. (2019) speculated that seamounts in the northern Ross Sea are the spawning grounds of D. mawsoni.The fishable habitats in the Ross Sea, including the spawning grounds in the deeper waters, were predicted to persistent under future environment conditions, highlighting the importance to protect these climate refugia from anthropogenic disturbance.

Figure 1
Figure 1 The predicted habitat suitability of Dissostichus mawsoni under current environment conditions.The available occurrence records were showed by pink points in the map.The colors, from blue to red, indicated the habitat suitability from low to high.The five lines indicated the five Antarctic Circumpolar Current fronts.NB, northern boundary; SAF, Subantarctic Front; PF, Polar Front; SACCF, Southern Antarctic Circumpolar Current Front; SB, southern boundary.Full-size  DOI: 10.7717/peerj.17131/fig-1

Figure 3 Figure 4
Figure 3 The ability of 10 individual models in predicting the habitat suitability of Dissostichus mawsoni, estimated by AUC and TSS values.TSS, True Skill Statistics; AUC, Area Under the receiver operating characteristic Curve.Data are expressed as mean ± standard error.Dashed lines indicated the cutoff values of TSS and AUC used to build ensemble model.Full-size  DOI: 10.7717/peerj.17131/fig-3

Table 1
Summary statistics of current environment conditions used to model the Dissostichus mawsoni distribution.

Table 2
Current environment conditions and the average values of climatic change in 2050s and 2100s under four Representative Concentration Pathways (RCPs) scenarios in the Southern Ocean.

Table 3
Change (%) of fishable distribution range of Dissostichus mawsoni under four Representative Concentration Pathways (RCPs) scenarios in 2050s and 2100s.